Clustering of extreme values: estimation and application

The extreme value theory (EVT) encompasses a set of methods that allow inferring about the risk inherent to various phenomena in the scope of economic, financial, actuarial, environmental, hydrological, climatic sciences, as well as various areas of engineering. In many situations the clustering effect of high values may have an impact on the risk of occurrence of extreme phenomena. For example, extreme temperatures that last over time and result in drought situations, the permanence of intense rains leading to floods, stock markets in successive falls and consequent catastrophic losses. The extremal index is a measure of EVT associated with the degree of clustering of extreme values. In many situations, and under certain conditions, it corresponds to the arithmetic inverse of the average size of high-value clusters. The estimation of the extremal index generally entails two sources of uncertainty: the level at which high observations are considered and the identification of clusters. There are several contributions in the literature on the estimation of the extremal index, including methodologies to overcome the aforementioned sources of uncertainty. In this work we will revisit several existing estimators, apply automatic choice methods, both for the threshold and for the clustering parameter, and compare the performance of the methods. We will end with an application to meteorological data.


Introduction
Climate change is at the order of the day, leading to a growing concern about the occurrence of phenomena such as extreme drought, floods, and large-scale forest fires. The pandemic situation caused by COVID-19 and the outbreak of armed conflicts as happened recently also have a strong impact on the global society and economy in which we live. Perhaps we have never seen an occurrence of extreme phenomena like today and a consequent demand for the use of appropriate tools to assess their impact, such as those provided by the extreme values theory (EVT).
The clustering effect of high values has a strong impact on the assessment of the risk associated with extreme phenomena and its main actor is the extremal index, denoted by . Indeed, describes and quantifies the clustering amount of the extreme values in many stationary time series. In Fig. 1 we can see the daily maximum temperatures collected at the climatologic station Abrantes in center of Portugal (in Celsius degrees). Clustering of extreme values is visible and thus the presence of extremal local dependence. Recently, other areas such as Dynamical Systems have also been applying the concept of extremal index (Moloney et al. 2019;Freitas et al. 2021, among others).
Let X X X = {X n } n≥1 be a stationary sequence of random variables (r.v.) with common marginal distribution function (d.f.) F. We say that X X X has extremal index ∈ [0, 1] if for each > 0 there exists a sequence of normalized levels u n , i.e., n(1 − F(u n )) → , as n → ∞ , such that P(M n ≤ u n ) → exp(− ) , where M n = max(X 1 , ..., X n ) . If = 1 then the tail behavior of X X X resembles an iid sequence, whenever < 1 leads to the occurrence of clusters of extreme values.
One extremal local dependence condition that is usually considered is the D(u n ) , which is basically a standard mixing condition that limits the long-range dependence at large values. It implies that any two exceedances of large u n , X i > u n and X j > u n , for sufficiently separated time points i and j are asymptotically independent (see, Leadbetter 1974).
If X X X satisfies D(u n ) , we have P(M n ≤ u n ) ≈ F n (u n ) , for large n and u n . Moreover, if there exist normalizing real constants a n > 0 and b n such that F n (a n x + b n ) → G(x) , then G is the d.f. of a generalized extreme value distribution (GEV) and P(M n ≤ a n x + b n ) → H(x) ≡ G (x) . If we consider {X * n } n≥1 an iid sequence with the same marginal d.f. F of X X X , the limiting GEV distribution of the corresponding result for M * n = max(X * 1 , ..., X * n ) is G(x) = H 1∕ (x) . In this context, plays a key role on the sample maxima distribution. The location and scale parameters of the GEV d.f. H and G are respectively related as follows: H = G + G ( − 1)∕ and H = G , where the shape parameter is the same in both G and H. Thus ignoring may lead to misspecified tail inferences: underestimation of quantiles of F if inference is based on H from sample block maxima or overestimation of H quantiles if inferences are based on marginal F from sample observations (see Beirlant et al. 2004). The extremal index also corresponds to the reciprocal of the mean cluster size in the point process of exceedance times of a large threshold u n , under a suitable mixing condition slightly stronger than D(u n ) (Hsing et al. 1988). Another interpretation of due to O'Brien (1987) is based on a conditional probability that quantifies to what extent extremes cluster together.
Many contributions on the extremal index estimation are addressed in the literature, based on different interpretations of . The major discussion focuses on strategies for the best choice of one or more auxiliary parameters involved in each method and respective stability of the estimates. Typically, the estimation of involves two sources of uncertainty: a threshold and some clustering parameter. Classical methods proposed in, e.g., Hsing (1991), Hsing (1993), Smith and Weissman (1994) and Weissman and Novak (1998) require both the choice of a threshold and a parameter associated with the clusters identification (runs parameter). More recently, threshold-dependent estimators have been proposed which are defined from the interexceedance times of a high threshold. This approach is based on the compound Poisson character of the point process of exceedances, namely, with an appropriate normalization and under a suitable and not restrictive mixing condition, the interexceedance times follow approximately an exponential mixture distribution with a point mass at zero and involving a parameter corresponding to . These include the intervals estimator of Ferro and Segers (2003), the K-gaps estimator of Süveges (2010), the censored and the truncated estimators of (respectively, Holěsovský and Fusek 2020;Holesovsky and Fusek 2022). The cycles estimator of Ferreira and Ferreira (2018) also requires the choice of a threshold and the validity of a local dependence condition describing the cluster behavior. The maxima estimators of Gomes (Gomes 1993;Ancona-Navarrete and Tawn 2000;Northrop 2015) are based on the comparison between the block maxima distribution obtained from the stationary sequence and the corresponding sequence of independent variables. We also include in this group the more recent estimator of Ferreira and Ferreira (2022) which is obtained from the bivariate block maxima distribution generated from the stationary sequence and an iid sequence of independent standard Fréchet variables. These estimators require the choice of a block size. In this paper we introduce a new version of the block maxima estimator in Ferreira and Ferreira (2022) and propose a bootstrap method to compute confidence intervals. In Sect. 2 we present a survey on the inferential methods that will be used. We are going to analyze their performance trough simulation in Sect. 3. Besides the pointwise estimation we also address interval estimation mainly based on bootstrap. In Sect. 4 we will illustrate the methods on a climatological dataset. We end with a discussion in Sect. 5.

Some characterization and estimation of
In this section we describe some probabilistic characteristics of the extremal index that inspired the methodology and the mathematical expression of the estimators to be presented. These will be analyzed on the next section through simulation and their performances will also be compared.
We start with classical runs estimator that is related to the O'Brien (1987) characterization, The empirical counterpart of (1) leads to the runs estimator where N u denotes the number of exceedances of threshold u in stationary sequence X X X and 1 {⋅} denotes an indicator function (Hsing 1993). In practice, clusters are identified by considering two different groups of exceedances of a threshold u as independent clusters if there are at least r − 1 consecutive observations below the threshold between them. We shall denote r the runs parameter. Hsing and McCormick (1991) establishes a similar result to that of O'Brien, under a local mixing condition denoted D (s) (u n ) which states that within a cluster, an exceedance of a high threshold u n is most likely to be followed by another exceedance within s − 1 consecutive observations. Condition D (s) (u n ) requires the validity of mixing condition D(u n ) that limits the long-range dependence at extreme levels by implying that any two exceedances of u n that are sufficiently separated in time are asymptotically independent. More precisely, X X X satisfies condition D(u n ) if for any integers 1 ≤ i 1 < ... < i q < j 1 < ... < j q ′ ≤ n for which j 1 − i q ≥ l , we have with n,l n → 0 , as n → ∞ , for some sequence l n = o(n) and l n → ∞.
Condition D (s) (u n ) will hold for X X X if D(u n ) also holds and there exists s > 0 integer, sequences r n and l n of integers such that r n → ∞ , n n,l n ∕r n → 0 , l n ∕r n → 0 and It is easily seen that once condition D (s) (u n ) holds, then D (s * ) (u n ) also holds for all s * ≥s.
Consider u n ( ) such that, for > 0, If D (s) (u n ) holds for the stationary sequence X X X , for some s > 0 and u n = u n ( ) for all > 0 , the extremal index of X X X exists if and only if for all > 0 ( Hsing and T., McCormick, W.P. 1991, Corollary 1.3).
Observe that the runs estimator (2) also corresponds to the empirical counterpart of (3) by considering the runs parameter r = s . When r = 2 we derive the Nandagopalan (1990) estimator which requires condition D (2) (u n ) to hold. Its formula is very easy to compute since it is simply the ratio between the number of upcrossings (or downcrossings) and the number of exceedances of u n . The approach in Ferreira and Ferreira (2018) is based on the Nadagopalan's estimator through an estimation procedure of the extremal index of an auxiliary stationary sequence satisfying D (2) (u n ) . More precisely, if X X X satisfies D (s) (u n ) , we take the so called cycles process {Z n = M (n−1)(s−1),n(s−1) } n≥1 for which D (2) (u n ) holds and estimate as the ratio between the number of upcrossings of threshold u within {Z 1 , ..., Z [n∕(s−1)] } , denoted U Z u and the number of exceedances N u of X X X (Ferreira and Ferreira 2018, Proposition 2.3), i.e., here denoted cycles estimator.
Consider the interexceedance times r.v. T(u n ) = min{j ≥ 1 ∶ X j+1 > u n |X 1 > u n } . Under a suitable mixing condition, we have that P(X 1 > u n )T(u n ) converges in distribution to a mixture distribution which is degenerated at zero with weight 1 − and has exponential law with mean value 1∕ with weight . Therefore, the extremal index expresses both the proportion of intra-cluster (within a cluster) times and inter-cluster (between clusters) times, and the expected value of the inter-cluster times under a convenient normalization. The intervals estimator of Ferro and Segers (2003) corresponds to a moment-based estimator derived from the limiting mixture distribution. It only requires the choice of a high threshold exempting the choice of a runs parameter. It will be denoted ̃( I) .
Although interexceedance times, T i = j i+1 − j i , i = 1, ..., N u , are not independent and a likelihood procedure assumes independence, this assumption may be disregarded under the validity of condition D (s) (u n ) for some s (Süveges 2007;Süveges 2010). On the other hand, the normalized intra-cluster times are theoretically zero in the limit but they are observed as positive values and thus will be assigned to the exponential part of the limiting mixture law.
By considering the new r.v. K-gap S (K) (u n ) = max(T(u n ) − K, 0) , the smallest times T i are set to zero, which improves the identification of intra-cluster times. Süveges (2010) shows that the maximum likelihood (ML) method can be applied to the limit mixture model in order to estimate , replacing interexceedance times by K-gaps S (K) and under condition D (s) (u n ) for s = K + 1 . In Süveges (2007), it was only addressed the case K = 1 . This method corresponds to the so called K-gaps estimator, here denoted ̃( K) . Choosing K involves some care, since if too large leads to more null times and consequent assignment to the degenerate part of the limiting mixture law, while too small K means fewer null times and hence an assignment tendency to the exponential component of the mixture law. The censored estimator introduced in Holěsovský and Fusek (2020) is similar to the K-gaps estimator but the smallest interexceedance times (less than some integer c) are censored. The choice of c is not so sensitive as the choice of K, however, the likelihood expression underlying the censored estimator does not allow an explicit formula for it, which makes its analysis and improvements difficult. A new approach is considered in Holesovsky and Fusek (2022), where the extremal index estimation is based on truncation. The small interexceeding times are truncated by a given parameter t, corresponding to K in the K-gaps estimator or c in the censored estimator. If condition D (s) (u n ) holds for X X X with s = t + 1 , under the same mixing condition considered to derive the limiting mixture law of P(X 1 > u n )T(u n ) , Holesovsky and Fusek (2022) prove that P(X 1 > u n )(T(u n ) − t)|T(u n ) > t converges in distribution to an exponential law with expected value 1∕ . Let T (1) ≤ ... ≤ T (N u −1) be the order statistics of T 1 , ..., with N t the number of times that are greater than some fixed positive t, and {S 1 , ..., S N t } = {T (N−N t ) − t, ..., T (N−1) − t} the set of exceedance times above the truncation value t. The ML method can be applied since truncated times are not affected by the sequence dependence (inter-cluster times are asymptotically independent (Ferro and Segers 2003) and intra-cluster times have not propensity to exceed t under local dependence condition D (t+1) (u n ) ). A simple ML estimator for corresponds to the arithmetic inverse of sample mean The derivation of the bias of (5) and of a penultimate approximation of the limiting distribution leads to an improved and bias corrected estimator, yielding the so-called truncated estimator with where ̃ is given in (5).
In applications, it can be difficult to check the validity of D (s) (u n ) condition. Various proposals have been presented, such as diagnostic plots of anti-D (s) (u n ) (Süveges 2007;Ferreira and Ferreira 2018), the information matrix test (Süveges 2010;Fukutome and Süveges 2014;Fukutome et al. 2019), or based on a stability check of the runs estimator (Cai 2019). However, the study of this issue is not closed and still awaits developments. The automation procedure in Fukutome and Süveges (2014); Fukutome et al. (2019) allows to select both s and u n . More precisely, considering the K-gaps estimator, it is based on misspecification tests through the information matrix test (IMT) presented in ( Süveges 2010). All combinations of pairs of thresholds and run parameters in plausible ranges are tested for misspecification of the model, and the pair (u, K) that generates the largest number of observations after declustering, within a list of pairs of small misspecification (IMT < 0.05) is selected, provided the number of exceedances is larger than 80. We will apply this automation procedure to select both threshold and runs parameter involved in each of the estimators: runs ̃( R) , cycles ̃( C) , intervals ̃( I) , K-gaps ̃( K) and truncated ̃( T) , in order to compare their performances through simulation and further on the analysis of real data. For a given selected pair (u, K), we will assume the validity of condition D (K+1) (u) . In the case of the intervals estimator, only the IMT threshold selection will be used since it solely depends on the threshold choice.
There are other estimation methods of entailing the choice of a single tuning parameter. Maxima methods as described in Gomes (Gomes 1993;Ancona-Navarrete and Tawn 2000;Northrop 2015) and more recently in Ferreira and Ferreira (2022) require the choice of a block length in order to generate a block maxima sequence from the original one. The two first references present methods that need to resample the original data to produce a sample of block maxima with approximate d.f. G. Northrop proposal in Northrop (2015) avoids this drawback by comparing the limiting GEV H of the maxima M n of X X X directly to the marginal d.f. F. More precisely, for large enough n, we have H ≈ F n , thus Y = −n log F(M n ) has exponential law with mean value 1∕ . The Northrop estimator is based on the ML approach, considering the sample of block maxima .., X n } and estimating the unknown d.f. F by the respective empirical d.f.. This corresponds to the disjoint blocks estimator. The sliding blocks version of Northrop estimator is based on a sample of overlapping block We will use the sliding blocks ̃( N) which compares favorably with the disjoint blocks (Northrop 2015). This will be denoted the Northrop estimator.
In Ferreira and Ferreira (2022) a new block maxima method was introduced to estimate the extremal index. It is based on a bivariate sequence {(Y n,1 = X * n , Y n,2 = (1∕2)X * n ∨ (1∕2)X n )} n generated from the original X X X which is assumed to have standard Fréchet marginals and by an i.i.d. sequence {X * n } n also having standard Fréchet marginals. Operator ∨ stands for "maximum". It is proved that the component-wise maxima of sequence {(Y n,1 , Y n,2 )} n . has a limiting bivariate extreme value copula C(u, v) = min(uv 1+ , v) with tail dependence coefficient leading to estimator

Algorithm:
The estimation procedure follows the steps below: Step 1. Take the marginal transformation − 1 where F X is an empirical d.f. of sample X 1 , ..., X n with dimension n, in order to have approximately standard Fréchet marginals.
Step 5. Repeat steps 2-4 a large number R of times, obtain estimates ̃ 1 , ...,̃ R and estimate in order to achieve robustness given the existence of arbitrariness in the generation of a random sample ( Step 2) in each estimate.
Step 4 is based on estimator (9 which requires an estimate of the tail dependence coefficient . We follow the proposal of Ferreira and Ferreira (2022) with where G j , j = 1, 2 , is an empirical distribution function of the GEV marginal G j of Z 1,j . For more details, see Ferreira 2022. In Ferreira andFerreira (2022) it was considered R = 10000 and block maxima taken on disjoint blocks. Here we consider the sliding blocks approach as in the Northrop estimator. Some prior simulations lead us to the proposal that R = 100 is reasonable for a robustness of the method. This approach will be denoted Ferreira estimator.
In the following we provide a catalog of some distributions and their associated extremal indices. Additional lists of models and respective extremal indices are also exposed in Gomes andNeves (2015, 2020).

Simulation study
In this section we analyze the estimation of the extremal index through simulation. We compare the performances of the runs (R) estimator in (2), the cycles (C) estimator in (4), the truncated (T) estimator in (6), the intervals (I) estimator of Ferro and Segers (2003), the K-gaps (K) estimator of Süveges (2010), along with the block maxima estimators of Northrop (N) in (7) and Ferreira (F) in (10). The first five estimators require the selection of a runs parameter and a threshold, except the intervals estimator only needing the choice of a threshold. To this end we apply the IMT method ( Süveges 2010;Fukutome and M.A. Süveges M. 2014;Fukutome et al. 2019) described in Sect. 2. The block maxima estimators require the settlement of a block length and no automation procedure is considered. In the study we take block lengths b = 10, 20, ..., 70 . The software( R Core Team 2020) was used and the R codes of estimators can be seen in https:// github. com/ msfer reira uminho/ msrf. The simulation study is based on the following models: a first order max autoregressive (MAR) process with standard Fréchet marginals and autoregressive parameter = 0.5 ( Davis and Resnick 1989) satisfying condition D (2) (u n ) (see, e.g.,Cai (2019)); a moving maxima (MMFrec) process with coefficients a j = 1∕5 , j ∈ {1, 2, 3, 7, 8} ( Deheuvels 1983) for which D (5) (u n ) holds (Ferreira and Ferreira 2018); a Markov chain (MCBEV) with standard Gumbel marginals and logistic joint distribution with dependence parameter = 0.5 (Smith 1992); an ARCH(1) process with Gaussian innovations, autoregressive parameter = 0.5 and variance parameter = 1.9 ⋅ 10 −5 (Embrechts et al. 1997); an AR(1) process with Cauchy marginals and auto-regressive parameter = −0.6 ( Chernick 1978) and a negatively correlated uniform AR(1) process with r = 2 ( Hsing and McCormick 1991), respectively denoted ARCau and ARUnif and both satisfying condition D (3) (u n ) . The extremal index values of the processes MAR, MMFrec, MCBEV, ARCH, ARCau and ARUnif are 0.5, 0.2, 0.328, 0.835, 0.64 and 0.75, respectively. We also consider a classical first order AR(1) process with Gaussian marginals and auto-regressive parameter 0.5, here denoted AR, which is almost independent and satisfies condition D (1) (u n ) having = 1 . As far as we know, there is no theoretical analysis on condition D (s) (u n ) for MCBEV model and we mention an empirical evaluation conducted in Ferreira and Ferreira (2018). Cai (2019) proved that D (s) (u n ) doesn't hold for a particular ARCH model.
We consider 1000 replicas of each model and compute the mean, the absolute bias (abias), the root mean squared error (rmse) and the standard deviation (sd). Bootstrap confidence intervals were obtained through percentile method technique on time series. Confidence intervals were obtained through bootstrap percentile method technique on time series with fixed block length 20 ( Kunsch 1989;Politis and Romano 1994), except in the case of the intervals estimator where the bootstrap proposal in Ferro and Segers (2003) was used. In the case of the likelihood estimators (K-gaps and Northrop block maxima), the confidence intervals are based on the asymptotic Normality. In Ferreira sliding blocks estimator we considered percentiles 2.5 and 97.5 applied on the auxiliary estimates ̃ j , j = 1, ..., R , produced by the method in Step 5 of the Algorithm in Sect. 2. We present the the proportion of intervals in the simulations that included the true value of (cov), the mean range width (rw) and the rate cov/rw corresponding to the proportion of coverage (cov) divided by the mean range width (rw). Some of the confidence intervals associated to Northrop sliding blocks failed to be computed because of the invalid standard error estimates which were accounted in column "NAs" of Table 2. The results for the runs, cycles, intervals, truncated and K-gaps estimators are presented in Table 1 and the estimates of Ferreira sliding blocks derived in (10) can be seen in Table 3. The rate cov/rw is plotted in Fig. 2 obtained for the runs, the cycles, the intervals, the truncated and the K-gaps estimators (left-top), for the Ferreira and Northrop sliding blocks estimators with block lengths b = 10, 20, ..., 70 (right-top) and the third bottom plot includes all estimators where for each of the sliding blocks estimators we address the best scenario (corresponding to the block length with the estimated largest rate, denoted by B) and the worst scenario (where the block length choice led to the smallest estimated rate, denoted by W). Analogous plots are represented in Figs. 3 and 4 relating to abias and rmse, respectively.
Cases closer to the border values of the extremal index domain were also considered, namely, MAR( = 0.9 ) with = 0.1 and MAR( = 0.1 ) with = 0.9 . In the almost independence Gaussian AR model with = 1 , besides parameter = 0.5 which induces strong dependence, we also analyze the AR weaker dependence models AR( = 0.1 ) and AR( = −0.1 ) where is still unit but they are even more close of independence. See Tables 4, 5 and 6.
The cycles (C), the runs (R) and K-gaps (K) estimators present a somewhat homogeneous performance across the various models, with the K-gaps estimator slightly better behaved. We recall that we are applying the IMT method in all estimators, except in the sliding blocks, which was developed under the K-gaps estimator. On Table 2 Simulation results obtained for the Northrop (2015) sliding blocks estimator: mean, absolute bias (abias), root mean squared error (rmse), standard deviation (sd), the coverage proportion of the true extremal index value within the estimated confidence intervals (cov), the intervals range width (rw), the ratio between "cov" and "rw" and the number of replicates where the confidence intervals couldn't be calculated ( the other hand, the intervals (I) and the truncated (T) estimators show some sensitivity to the different models: they perform very well in the MMFrec, MAR, MCBEV and ARCH, but their biases and rmse are high in the ARUnif model. We can also see that the choice of the block-maxima size to be used in sliding blocks estimators is important as their behavior improve in the best block choices. Indeed, both Ferreira and Northrop sliding blocks estimators present the best overall performances along the different models under the best block choices. Both Northrop and Ferreira sliding blocks have to deal with a trade-off between bias (larger for small block-maxima lengths) and variance (increases with the block-maxima length). This can be observed in Tables 2 and 3. In practice one can plot estimates for several block-maxima sizes b and choose the smallest b above which estimates are approximately constant (Northrop 2015).
In the coverage rate (Fig. 2), it is observed that the K-gaps estimator presents an overall better performance. Yet, the coverage percentages ("cov") in Table 1 are not as large as would be desirable. This is particularly notorious in cases where estimates present larger rmse like AR model. In order to analyze the effect of bootstrap block-size choice on CI coverage, we have also applied the automatic block-size choice method developed in Politis and White (2004) and (Patton et al. 2009), available in R package blocklength Stashevsky 2022. We recall that the bootstrap CI requiring block-size selection was considered for the runs, cycles and truncated estimators. The results of the automatic block-size choice method are in Table 7 where we observe improvements for the truncated estimator, except in the ARCH model. On the other hand, the ARCH model benefits from the automatic method within the runs and cycles estimators. The  Table 3 Simulation results obtained for the ( Ferreira and Ferreira 2022) sliding blocks estimator: mean, absolute bias (abias), root mean squared error (rmse), standard deviation (sd), the coverage proportion of the true extremal index value within the estimated confidence intervals (cov), the intervals range width (rw), the ratio between "cov" and "rw" and the number of replicates where the confidence intervals couldn't be calculated ( MAR model with = 0.1 also gains with the automated method. The Ferreira sliding blocks estimator presents both the highest coverage proportions (cov) and range widths (rw), leading to low coverage rates ("cov/rw"), as can be seen in Table 3 and Fig. 2. The compromise between both high coverage ("cov") and coverage rate ("cov/ rw") points out the Northrop sliding blocks estimator as the best choice, particularly for larger block-maxima sizes. The AR model with stronger dependence (i.e., dependence parameter = 0.5 ) presents the worst performance in all estimators (see, e.g., Figs. 2, 3 and 4). Both sliding blocks estimators seem to be the most promising in this model, which has an extremal index equal to one and, therefore, a boundary value in the support of . However, the same is not true with weaker dependent AR models (i.e., with dependence parameter = ±0.1 ), where the absolute bias and rmse are lower, particularly within intervals (I), truncated (T) and Northrop (N) estimators (Tables 4, 5 and 6). The bootstrap confidence intervals of runs (R) and cycles (C) exhibit very small coverage percentages as well as the ones of K-gaps estimator based on ML, for all considered AR models. In the challenging cases of with values close to the domain boundaries 0 and 1 through model MAR with autoregressive parameters = 0.1 and = 0.9 , corresponding to = 0.9 and = 0.1 , respectively, we observe that estimators of runs (R), cycles (C), K-gaps and Ferreira sliding blocks (F) tend to behave better for = 0.1 than = 0.9 while in the truncated (T) estimator the conclusion is opposite. The intervals (I) and Northrop sliding blocks (N) present an overall better performance in both situations. Table 4 Simulation results of models MAR with = 0.1 and = 0.9 , AR( = 0.1 ) and AR( = −0.1 ), obtained for the runs (R), cycles (C), intervals (I), truncated (T) and K-gaps (K): mean, absolute bias (abias), root mean squared error (rmse), standard deviation (sd), the coverage proportion of the true extremal index value within the estimated confidence intervals (cov), the intervals range width (rw) and the ratio between "cov" and "rw"

Mean
Abias  Table 5 Simulation results of models MAR with = 0.1 and = 0.9 , AR( = 0.1 ) and AR( = −0.1 ), obtained for the ( Northrop 2015) sliding blocks estimator: mean, absolute bias (abias), root mean squared error (rmse), standard deviation (sd), the coverage proportion of the true extremal index value within the estimated confidence intervals (cov), the intervals range width (rw), the ratio between "cov" and "rw" and the number of replicates where the confidence intervals couldn't be calculated ( Table 6 Simulation results of models MAR with = 0.1 and = 0.9 , AR( = 0.1 ) and AR( = −0.1 ), obtained for the Ferreira and Ferreira (2022) sliding blocks estimator: mean, absolute bias (abias), root mean squared error (rmse), standard deviation (sd), the coverage proportion of the true extremal index value within the estimated confidence intervals (cov), the intervals range width (rw), the ratio between "cov" and "rw" and the number of replicates where the confidence intervals couldn't be calculated (

Application
We consider the daily maximum temperatures in the months of July and August, from 2001 to 2021, collected at the climatologic station Abrantes in center of Portugal (in Celsius degrees) available at "SNIRH: Sistema Nacional de Informação dos Recursos Hídricos". 1 The data is plotted in Fig. 5, in which successive high temperatures are seen, which are typical in this inland city. The results in Table 8 are obtained through the application of IMT procedure, leading to threshold 37.7 and K = 2 . Thus we assume D (3) (u n ) dependence condition. Figure 6 represents the estimates of ̃( R) , ̃( I) , ̃( C) . ̃( T) and ̃( K) , for thresholds corresponding to quantiles from 0.8 to 0.99. The estimators give quite close results to each other, except the sliding blocks ̃( F) with larger values (Fig. 6, left). However, most of ̃( F) estimates are within the Northrop sliding blocks confidence bands (Fig. 6, right), which we use as reference according to the simulation study findings in Sect. 3. Our guess is that possible values for range between 0.4 and 0.55.

Discussion
The extremal index is a very important measure in inferring extreme values of time series. In addition to affecting the limiting law behavior of the maximum, it is also associated with a clustering effect of exceedances of high values that can have Table 7 Coverage proportion of the true extremal index value within the estimated bootstrap confidence intervals (cov), the intervals range width (rw) and the ratio between "cov" and "rw", considering automated choice of block length in Politis and White (2004) and Patton et al. (2009)  harmful consequences. The estimation of the extremal index is a long-standing subject in extreme value theory, and it still attracts the attention of researchers. Here we present a study that covers several estimators, from the classic runs estimator to the most recent methodologies in block maxima, such as ( Ferreira and Ferreira (2022)) work resorting the theory of bivariate extremes. This new approach is still at an early Fig. 2 The ratio between the coverage proportion of the true extremal index value within the estimated intervals and the respective range width ("cover/range") obtained for: (left-top) the runs (R), the cycles (C), the intervals (I), the truncated (T) and the K-gaps (K); (right-top) the sliding blocks of ( Ferreira and Ferreira (2022)) (black) and of ( Northrop (2015), ) (grey), for block lengths b = 10, 20, ..., 70 ; (bottom) the runs (R), the cycles (C), the intervals (I), the truncated (T) and the K-gaps (K), the sliding blocks of ( Ferreira and Ferreira (2022)) (black) and of ( Northrop (2015)) (grey), for block lengths corresponding to the largest "cover/range", denoted by B (from best) and with the smallest "cover/range", denoted by W (from worst) stage, and there is room to improve both the estimation methodology that requires the generation of auxiliary samples and the analysis of the asymptotic behavior of the estimator and obtaining confidence intervals. Finding the asymptotic variance of the extremal index estimators is a great challenge given the context of serial dependence. Thus the resampling methodologies appear as good alternatives. However the Fig. 3 The absolute bias (abias) obtained for: (left-top) the runs (R), the cycles (C), the intervals (I), the truncated (T) and the K-gaps (K); (right-top) the sliding blocks of (Ferreira and Ferreira (2022)) (black) and of Northrop (2015) (grey), for block lengths b = 10, 20, ..., 70 ; (bottom) the runs (R), the cycles (C), the intervals (I), the truncated (T) and the K-gaps (K), the sliding blocks of Ferreira and Ferreira (2022) (black) and of Northrop (2015) (grey), for block lengths corresponding to the smallest "abias", denoted by B (from best) and with the largest "abias", denoted by W (from worst) choice of bootstrap block length is a critical point of estimation. A careful analysis is still needed to ameliorate in practice the coverage probabilities which proved to be poor in some cases. The great advanced in computing capabilities that we are currently witnessing opens up good prospects for the implementation and development of techniques such as bootstrap or jackknife. The area of machine learning is also a whole new horizon waiting to be explored.  (2015) (grey), for block lengths corresponding to the smallest "rmse", denoted by B (from best) and with the largest "rmse", denoted by W (from worst)